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I. INTRODUCTION 



We consider the eigenfrequency spectrum of isotropic and homogeneous elastic bodies of arbitrary shape in two 
dimensions. The governing linear equations, the Navier-Cauchy equations, are separable only for a very small set of 
geometries such as spherical bodies or infinitely long cylindrical wave guides. Solutions to the vast majorities of shapes 
can be obtained only with the help of numerical techniques such as finite element or boundary integral methods 0,0] . 
Purely numerical approaches are, however, severely limited by computer resources and often restricted to the low 
frequency regime with wavelengths only one or two orders of magnitudes smaller than the typical size of the system. 
In the high frequency limit statistical methods such as statistical energy analysis Q or random matrix theory £§] have 
proved valuable. While the former yields information about mean response signals neglecting interference effects, the 
latter provide answers regarding the universal part of the fluctuations in the signal not taking into account system 
dependent effects. An alternative approach providing more detailed information in the mid to high frequency regime 
is obtained using asymptotic methods. 

Similar to geometric optics, asymptotic methods connect wave propagation to an underlying ray dynamics in the 
■ small wave lengths limit. In elasticity, typical boundary conditions such as free boundaries lead to mode coupling and 
' ray splitting at boundaries which complicate a ray analysis. The bulk of the asymptotic analysis in the elastodynamics 
literature has focused on interface or scattering problems [f| in which the number of reflections for a typical path is 
small. A ray-treatment of interior modes of elastic bodies of finite size has received much less attention so far; here, 
geometric rays have an infinite number of reflections with the boundaries undergoing ray-splitting at every impact 
which leads to summation and ordering problems when expressing operators such as the Green function in a short 
wavelength approximation. 

Such issues have been addressed in the context of the Helmholtz equation in finite domains and more generally in 
quantum mechanics. Especially, the connection between the solution of these scalar wave equations and the dynamical 
joperties of a related ray or classical dynamics have been treated in much detail in the context of quantum chaos 
■ A powerful tool connecting the spectrum of a quantum system with an underlying classical dynamics are 
trace formulae as for example introduced by Gutzwiller in quantum mechanics |(| , expressing the trace of the Green 
^ ' function in terms of the periodic orbits of the classical system. 

A trace formula for the interior problem in elasticity has been presented first in the result was, however, obtained 
by way of comparison with the scalar Helmholtz equation and not derived from the governing equations. To the best 
of our knowledge, such a derivation is still lacking, which is desirable not only from a point of principle, but is essential 
to obtain corrections beyond the leading large wavenumber asymptotics. A powerful technique for deriving such trace 
formulae from the underlying wave equations is to derive asymptotic expressions for boundary integral operators in 
terms of semiclassical transfer operators, a method pioneered by Bogomolny llfi , see also [Ullialia. The general 
form of such a transfer operator for elastic problems has been postulated in |l4| and verified for the special case of 
a circular waveguide. A derivation of the transfer operator for the biharmonic equation describing the out-of-plane 
vibrations of plates has been obtained in |15| incorporating the coupling of flexural and boundary modes. In the short 
wave length limit, the wave equation reduces again to a scalar problem with modified boundary conditions due to the 
exponential damping of surface waves away from the boundary. 

In the following we derive the Bogomolny transfer operator and form there periodic orbit trace formulae for two 
dimensional elasticity starting from first principles. 
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II. THE TRANSFER OPERATOR 
A. Fundamental equations 

We consider isotropic and homogeneous elastic bodies described in the frequency domain by the Navier-Cauchy 
equation |16( 

/iAu+(A + ^)V(V-u) + pw 2 u = 0, (1) 

where u(r) is the displacement field, A, fi are the material dependent Lame coefficients and p is the density which we 
assume to a constant. We will consider free boundary conditions here, that is no forces act normal to the boundary; 
this can be expressed in terms of the traction t(u), that is, 

t(u) = n • er(u) = (2) 

where n is the normal at r on the boundary C of the elastic body; the stress tensor <r(u) is given as 

<x(u) = A(V • u) 1 +/x(V® u + V) . (3) 

We make the standard Helmholtz decomposition of the displacement field u, that is, 

u = Up + u s with Up = V$, u s = V x * ; (4) 

the elastic potentials $ for the pressure (or longitudinal) and for the shear (or transversal) wave component solve 
Helmholtz's equation 

(A + fcjp$ = 

(A + fcf = 1 > 

with wave numbers k p and k s , respectively. One finds the dispersion relation fc PjS = uj/c p _ s with wave velocities 



/A + 2^ fjr 
V p Vp 

In the following, we shall restrict ourselves to two-dimensional problems, that is r, u(r) £ M 2 and we set \& = 
(0, 0, \E , ) t . The resulting differential equations describe in-plane deformations in plates or wave propagation in bodies 
with fixed shape in the xy plane extending to ±oo along z. 

B. Boundary integral equations 

1. General set-up 

In what follows, we will adapt the method outlined in [Tol Il5| to the Navier-Cauchy Eqn. Q). We first rewrite 
the boundary conditions J2J in terms of boundary integral equations and then consider the Fourier coefficients of the 
boundary integral functions. We start by introducing the elastic potentials in the form 



$(r) = J^G(r,a;k p )g(a)da; (7) 

*(r) = J G(r,a;k s )h(a)da. (8) 

where g and h are yet unknown single layer distributions on the boundary and a € [0, Lc] parameterises the boundary 
of length Lc, that is, r(a) S C; furthermore, G(r,r';fc) is a Green function solving the inhomogeneous Helmholtz 
equation 



(A + fc 2 )G(r,r';fc) = 5(r - r') 
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FIG. 1: Coordinates on the boundary: a) position representation with path of length L from r'(a) to r(/3) (here for a initial 
pressure wave); b) momentum representation with shear and pressure path starting with tangential momentum p' and ending 
with momentum p on the boundary. 



The integrals converge for r inside C and non-singular layer distributions g and h, and the ansatz J7J), J5J) thus solves 
the Helmholtz equation in the interior. A convenient choice for G(r, r';fc) is the free Green function which in 2 
dimensions takes the form 

G(r,r',k) = ^H^(k\r-r'\) (9) 

where Hq is the 0-th order Hankel function. 

In a next step, it is useful to rewrite the boundary condition J5J in terms of the elastic potentials. Defining n and 
t as the (outward) normal and tangent vectors at the boundary point v(/3) S C as indicated in Fig.^i, one obtains 



n • cr = -\hk 2 $ + 2/i 



d 2 , , d 2 , „ d 2 , f 1 _ d 2 \ 



n <I> + t $ -I- n ^ + t 

<9n 2 dndt dndt 2 V <9i 2 dn 2 I 



= 0, (10) 



where we used A$ = — k^Q valid in the interior; note that all partial derivatives are understood as being taken in 
the interior (after a suitable continuation of the local coordinates system into the interior) and then taking the limit 
r -> r(p) e C. 

We thus need to determine derivatives of the form 



d nn / G(J3,a)f{a)da; d nt / G(fi,a)f(a) da; 8 tt / G{fi, a)f(a) da (11) 
J C J c «c 

with G((3, a) = G(r(/3), r'(a), fc), / stands for g or h, respectively, and the derivatives are always taken with respect to 
the first variable r((3) from the interior. Note that taking the limit r — ► r((3) and differentiating are non-commuting 
operations due to the logarithmic singularity of the Green function for [3 — > a. 

For a short wavelength analysis, we distinguish between long segments with fc|r(/3) — r'(a)| 3> 1 and short con- 
tributions with fc|r(/3) — r'(a)| = 0(1); for the former, one can employ the asymptotic form of the Green function 



G(r,r',fc)~l J 2 e *(fc|r-rW4) fc | r _ r '|^ cx , j (12) 
4i y 7rfc|r — r'| 

whereas the logarithmic singularity G(r,r',k) ~ ^ln(fc|r — r'|) for k\r — r' | — ^ calls for a separate treatment for 
short length contributions. We note in particular, that one obtains from i|12|) in leading order, 

d n G(fl,a)~iqG{fl,a); d t G(/3,a) ~ ipG(p,a), (13) 

and likewise for the second order derivatives. Here, q{(3,a) — fccos^o and p((3,a) = fcsin^o are the normal and 
tangential component of the wave vector k = k (r(/3) — r'(a))/|r(/3) — r'(a)| at the boundary point j3, see Fig. 
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2. Asymptotic form of the boundary integral kernel in momentum representation 
Following |l5p. we split the boundary integral into two parts, that is, 



da = da + da 

JC/A J A 

where A refers to a small interval around a — (3 scaling as A ~ fc~ 1+e with < e < 1. 

We deal with the short length contributions first. Due to the scaling chosen for the interval A, we can neglect 
curvature contributions in the large k limit and write in leading order in 1/fc 

-i pkA/2 -i />oo 

G(r(/3),r(a),k)f(a)da ~ -- / G(0, x/k, k)f(x/k)dx ~ -- / G{0,x/k,k)f(x/k)dx (14) 

A K J-kA/2 * J-oo 

thus integrating along a straight line in direction of t(/3) centred at r(/3). It is now convenient to express the free 
Green function in integral representation, which in two dimensions leads to 

f dr> 2 e i P'< r - r ') 

"''■^'--a/s v-^+fcr (15) 

Aligning the x-axis with the tangential direction t(/3) as in Ijl4|l and integrating out the p y component, one obtains 

f dm t \ e i( i\y-y'\ 



with q = y/k 2 — p 2 . Note that lim^O- d y G(r, r') = ^S(x — x') revealing the singular behaviour of the Green function 
in this limit. 

Next, we express the single layer distributions on the boundary in its Fourier components, that is, 

f(a) = J dpf p e* pa (17) 

where we treat p to leading order as a continuous variable neglecting the discreteness of p = 2irj / Lc, jeN due to the 
finite length of the boundary C. From (|16|l together with (|14f) . we obtain the short length contributions in the form 

G(8,a)f(a)da~ - lim \ dpdp' e ipfi - / — e^ p - p ^ x ' k f p , = / dp^- e 1 *" 3 . (18) 

A k v^°- J 2 *9 J-oo 27r ' J 2 *9 

We proceed as above for the partial derivatives (|ll|l by identifying the normal and tangential direction with the y 
and x axis, respectively. Note again, that the derivatives d/dy need to be taken before completing the limit y —> 0_ 
from below. One obtains 

dnnj G(P,a)f(a)da = i J dj>/ p |e^; (19) 

d nt J^G(f3,a)f(a)da = -i J dpf p ^; (20) 

duj G{p,a)f(a)da = i Jdpf p ^-e^. (21) 

Turning to the contributions from long trajectories, we again introduce the Green function on the boundary in 
terms of its Fourier components 

G(/3, a, fc p/s ) = J dpdp' G^e^-^ (22) 

and write 

da G(f3, a)f(a) = [ dpdp' G pp > f p >e ipl3 . 

CI A J 
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Here, differentiation can be pulled under the integral sign and by employing the asymptotic form 11 2[) together with 
a stationary phase approximation, one obtains in leading order 

dnGppi — ipGppi; dtG P pi — iqG pp i 

and likewise 

dnnGppi — p G pp ' , dtnGpp' — qpGpp' , dttGpp' — q Gpp 1 . 

Writing the boundary conditions @ in terms of the Fourier components g p , f p and G^f for both short and long 
contributions, one obtains the set of equations 



(m + MD) X = with X p = f | p j (23) 



where 



and 



(M0) -' " 2 1 -2w %(p 2 ~ & J PP ' 1 PP ' ~ 2 { f s (P 2 - <&) J SpP ' (24) 



The eigenfrequency condition for finite elastic bodies in 2 dimensions can thus be cast into the form 

det(I - f (uj)) = with T = -Mg 1 MD (26) 



and 



«" 4det(M )l v ^(A^ + 2^2) ^(p 2 - g*)(Afc* + 2/^) + 4 M V J ' ^ lj 



as well as 



det(Mo) = -- 



^(p 2 -a(A^ + 2^)-4 M V 
9s 9p 



(28) 



The operator T is the short wavelength approximation of a wave propagator acting on boundary functions in Fourier 
or momentum representation; it has the general form of a quantum Poincare map [Tfl Il3j |. here written for the 
clastodynamic case including mode conversion. The matrix elements T pp / describe the evolution of pressure and shear 
waves along 'ray' - trajectories starting on the boundary with tangential momentum p' and hitting the boundary with 
tangential momentum p\ note that the rays corresponding to two different modes will in general start and end at 
different points on the boundary, see Fig. \I]p. The q p / s component is the part of the wave vector k p / s normal to the 
interface and we may set 

Pp/s = kp/s shiflp/s; q p/s = k p/s cos6» p/s . 

The tangential momentum p at the end points is the same for both polarisations before and after impact with the 
boundary and we obtain directly Snell's law 

p = p p = k p sin 6 p = fc s sin 9 S — p s . (29) 

Using k = k s /k p = c p /c s and identities like 

Xkl + 2/ic? 2 = (A + 2/i)fc 2 cos 29 s ; p 2 - q 2 s = -k 2 s cos 26 a , 

we may write the pre-factor matrix in the form 

A=-M- 1 M= (j^> A £ V (30) 
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with 



sin 29 s sin 26> p — n 2 cos 2 29 s 
pp ss = sin20 s sin2(? p + K 2 cos 2 2# s ( } 



2 2 sin 20 s cos 2# s 
p " K sin29 s sm29 p + k 2 cos 2 26 s 
2 sin 2# p cos 29 s 
~ sin 2<9 S sin 2(9 p + k 2 cos 2 26> s 



The matrix elements of A are up to a similarity transformation equivalent to the standard conversion factors for plane 
shear or pressure waves at impact with a plain interface and free boundary conditions |16j . Note that we follow here 
the convention used throughout the paper; for example, A sp denotes the conversion amplitude between an incoming 
p - wave and an outgoing s - wave. 

Next, we express the transition matrix A in a slightly different form using the transformation 

a = K-AK with K=(<*/f' 4 ) <*> 

which leads to a unitary matrix a. The relations a 2 + a 2 p = 1 = a ps + a 2 s reflect conservation of wave energy normal 
to the surface in the presence of mode conversion [l6j| . 



3. Asymptotic form of the boundary integral kernel in position representation 

It is often convenient to work with the boundary integral kernel in position representation; the inverse Fourier 
transformation of the operator T pp < = A p D P p< again taken in stationary phase approximation and employing the 
asymptotic form of the free Green function l|12|). yields 



1 a a fa \ ( v / fcpV fc ? L 



T(/3, a) = - 1 == cos 9 A(/3, a) ( V "v" " ) . (33) 



The stationary phase condition picks out contributions from shear and pressure waves travelling from a to /3 along 
rays of length L intersecting the boundary at (3 with a common angle (9 , see Fig.^,. In contrast to the momentum - 
representation considered earlier, rays leaving the end point (3 can do so along three different directions with angles 
9q, 9 P and 9 S - A p - polarised wave, for example, may emerge from (3 at an angle #o or 9 P depending on whether 
the corresponding incoming wave was a p or s wave. We thus set 9 p = 9q in A pp and A sp and 9 S = 9q in ^4 ps and 
A ss in Eqn. with 9 P ,9 S given by Snell's law l|29(l : note that this implies for example that A pp ^ A ss in general. 
Rewriting the operator (|33|l in terms of the (now in general non-unitary) transition matrix a, one obtains 



\j2niL V V K cos cos 9 S a sp cos 9 a ss J \ \fk~ s e" 



with 



T((5 a) = __ | cos6» a P p y/cos9 Q cos 9 p /k a ps \ _ f ^/k p e" ■ " k L I (34) 



(35) 



sin 29 s sin 29 Q - k 2 cos 2 29 s 2 Jam 29 sin 26» p cos 29 

— —K 



sin 29 s sin 29 + k 2 cos 2 29 s ' p sin 26» sin 29 p + k 2 cos 2 29 q 

2 Vsin 20 o sin 26* s cos 29 s sin 26>o sin 29 p — k 2 cos 2 20 o 



flsp K sin 26» s sin 29 + k 2 cos 2 26» s ' " ss sin 26» sin 29 p + k 2 cos 2 26> ' 

For hyperbolic shapes, that is, for boundaries only admitting isolated periodic geometric rays (including mode 
conversion at the boundary), standard arguments lead to a description of the traces of the operator T in terms of 
periodic ray trajectories [T(|. One obtains 

(n) 

TrT™ = Aje lS ^^l 2 (36) 

3 

where the sum is over all periodic ray trajectories having n reflections at the boundary with position and polarisations 
[(o^ , ) , . . . , {ct J n ,l J n )] where l\ = p or s is the polarisation of the ith segment of the periodic ray j leaving the boundary 
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at the point a,, i = 1, . . . , n. Furthermore, one has 



A — A geo 



i=l 



(37) 



taken along a periodic orbit; here Sj is the action of classical mechanics and the amplitude Aj separates into a 
geometric part Aj £ ° containing information about the spreading of nearby trajectories and a mode conversion loss 
factor. The traces T n contain all the information about the spectrum and may be used to construct the density of 
states or express the spectral determinant (|26f) . 
The operator 

can be written in a form more familiar from semiclassical quantum mechanics. We note that the 
cosine terms in the amplitudes relate to ray angles before and after hitting the boundary at [3; each contribution to 
the periodic orbit formula (|3fci|) thus contain products of cosine terms along the periodic orbit. Following an argument 
by Boasman ^3 developed in the scalar case, we consider 



d 2 L(P,a) 



dad (3 



cos 6 a p cos 9 



0a 



(38) 



with angles 9p a = 9q taken at (3 and 9 a p taken at a, respectively, (see Fig. ^,) 
IJ3D and f defined as 



The traces of the operators T as in 



T(/3,a) 



_i 



3 2 Lfl 



dad/3 



a(/3,a) 



'fcp e 



ik.Lr- 



(39) 



are thus equivalent to leading order. That is, when writing the traces as sum over periodic rays as in (13611 . the cosine 
terms coincide after multiplication along a periodic orbit. Similarly, the extra k ±:l / 2 terms in the off-diagonal terms 
in (|34|) cancel after one period. This confirms the form of the operator as postulated in 0] from which the trace 
formula suggested by Couchman et al can be derived by standard means as indicated earlier. 



III. CONCLUSION 



We have derived an asymptotic form of the boundary integral kernel in 2d elastodynamics from which periodic 
orbit trace formulae can be deduced using stationary phase arguments. It is expected that a 3d version of the 
asymptotic operator can be written in the form (|39|l using local coordinates where the tangential direction lies in 
the plane spanned by the vector r — r' and the normal at the boundary point r. In deriving the 3d version of the 
operator l|39|l one is naturally lead to a momentum representation in terms of spherical coordinates; the technical 
difficulties are not expected to exceed these of the 3d quantum case as discussed in ^(j and ES ■ 
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